# konrad.smolinski@gmail.com
# date: 26/11/2010
# last update: 1/12/2010
# info:
#
#
# ------------------------------------------------------------------
# Call for libraries:
rm(list=ls())
source("./rcLibSrc.R")
#Set Up
K <- 2 
# set up:
a <- c(0,1)
b <- c(0,1)
d <- c(2,3^(-1/2),3^(-1/2),3^(-1/2))
zval <- c(-0.5,0.5)
R <- length(zval)
xthres <- c(0)
xval <- c(0,1)

# integration 
acc <- 40
# Validate :
dgp <- rcDgp(a,b,d,xthres,xval,zval)
dgp

# this part is done once:
probIdx <- rcLeaveOutToProbs(K)

probMat <- matrix(0,ncol=R,nrow=length(probIdx))

for(r in 1:R){
	for(i in 1:length(probIdx)){
		 probMat[i,r] <- sum(dgp[[r]][probIdx[[i]]]) 
	}
}
rhs <- apply(probMat,1,max)



# integrand over Base Components
# this part depends on initial parameters
loIdx <- rcLeaveOutIdx(K)
lhs <- rep(0,length(loIdx))

# parameters enter here:
intVec <- rcBaseIntegrate(a=c(0,1.05),b=c(0,1),xval=xval)
for(i in 1:length(loIdx))	lhs[i] <- sum(intVec[loIdx[[i]]])
lhs

lhs > rhs

all(lhs > rhs)

# GRID
a0seq <- seq(-0.2,0.5,len=61)
a1seq <- seq(0.1,1.6,len=61)

aGrid <- expand.grid(a0seq,a1seq)

b2vec <- c(0.4,0.8,1,1.2,1.6,1.8,2.2)
res <- matrix(FALSE,nrow=nrow(aGrid),ncol=length(b2vec))

qs <- smolyak.quad(2,40)
set <- list()
polyAlpha <- list()

for(l in 1:length(b2vec)){
	cat('\n')
	cat(paste('b=',b2vec[l],sep=''))
	for(g in 1:nrow(aGrid)){
		lhs <- rep(0,length(loIdx))
		intVec <- rcBaseIntegrate(a=c(aGrid[g,][[1]],aGrid[g,][[2]]),b=c(0,b2vec[l]),xval=xval,quadPt=qs)

		for(i in 1:length(loIdx))	lhs[i] <- sum(intVec[loIdx[[i]]])
	
		 res[g,l] <- all(lhs > rhs)
		cat(".")
	}
	
	set[[l]] <- aGrid[res[,l],]
#	polyAlpha[[l]] <- ahull(set[[l]],alpha=0.2)
}


save(set,dgp,b2vec,a0seq,a1seq,a,b,d,zval,xval,xthres,file='../Data/rcEx1-1-Data.RData')
load(file='../Data/rcEx1-1-Data.RData')
# plot
polyAlpha <- list()

pdf("../Graphs/rxEx1-1.pdf",width=10,height=6)
#postscript("../Graphs/rxEx1-1.eps",width=10,height=6)
polyAlpha <- list()
par(mfrow=c(2,3))
for(i in c(1,3,4,5,6,7)){
	print(i)
	polyAlpha[[i]] <- ahull(set[[i]],alpha=0.1)
	plot(polyAlpha[[i]],main=paste("b[2]=",b2vec[i],sep="")
	,xlab=expression(a[1])
	,ylab=expression(a[2])
	,wlines='vor',col=c(1,2,3,3,0,2))
	
	abline(v=0,col="red",lty=1)
	abline(h=1,col="red",lty=1)
	points(0,1,pch=19,cex=1.2,col=2)
}


dev.off()